Load libraries

#library(GmAMisc) 
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readxl)
library(ggplot2)
#library(dplyr)
library(ggthemes)
library(RColorBrewer)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(ordinal) 
## 
## Attaching package: 'ordinal'
## The following object is masked from 'package:plotly':
## 
##     slice
## The following object is masked from 'package:dplyr':
## 
##     slice
library(emmeans)
library(multcompView)
library(broom.mixed) 
## Registered S3 method overwritten by 'broom.mixed':
##   method      from 
##   tidy.gamlss broom
library(forcats)
#library(plotly)
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:car':
## 
##     logit
library(egg)
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Bootstrapping functions

# functions for bootstrapping 95% confidence intervals around the mean

# get mean and 95% CI of values x via bootstrapping
boot_ci <- function(x, perms=5000, bca=F) {
  get_mean <- function(x, d) {
    return(mean(x[d]))
  } 
  x <- as.vector(na.omit(x))
  mean <- mean(x)
  if(bca){
    boot <- boot.ci(boot(data=x, 
                         statistic=get_mean, 
                         R=perms, 
                         parallel = "multicore", 
                         ncpus = 4), 
                    type="bca")
    low <- boot$bca[1,4]
    high <- boot$bca[1,5] 
  }else{
    boot <- boot.ci(boot(data=x, 
                         statistic=get_mean, 
                         R=perms, 
                         parallel = "multicore", 
                         ncpus = 4), 
                    type="perc")
    low <- boot$perc[1,4]
    high <- boot$perc[1,5] 
  }
  c(low=low,mean=mean,high=high, N=round(length(x)))
}


# get mean and 95% CI via bootstrapping of values y within grouping variable x
boot_ci2 <- function(d=d, y=d$y, x=d$x, perms=5000, bca=F){
  df <- data.frame(effect=unique(x))
  df$low <- NA
  df$mean <- NA
  df$high <- NA
  df$n.obs <- NA
  for (i in 1:nrow(df)) {
    ys <- y[which(x==df$effect[i])] #pulls out all the e.g.location prefs
    if (length(ys)>1 & var(ys)>0 ){
      b <- boot_ci(y[which(x==df$effect[i])], perms=perms, bca=bca) #resamples with replacement the mean and ci for e.g. AJs pref for location, gives back the low mean and high
      df$low[i] <- b[1]
      df$mean[i] <- b[2]
      df$high[i] <- b[3]
      df$n.obs[i] <- b[4]
    }else{
      df$low[i] <- min(ys)
      df$mean[i] <- mean(ys)
      df$high[i] <- max(ys)
      df$n.obs[i] <- length(ys)
    }
  }
  df
}
# y = response score, x= response within group (e.g. bat spp)

Import data:

#use relevant working directory


# gleaner_universe_r_7_19_19 <- read_excel("~/Dropbox/GLEANER UNIVERSE STUDY/data sheets/gleaner_universe_r_7.19.19.xlsx", 
#     na = "NA")





g_u <- read.csv("https://raw.githubusercontent.com/maydixon/katydids/master/gleaner_universe_r_08.15.22_github.csv", sep=",", header=TRUE, stringsAsFactors = F )



str(g_u)
## 'data.frame':    1933 obs. of  28 variables:
##  $ Species                     : chr  "Trin_nice" "Lonc_auri" "Micr_hirs" "Phyl_hast" ...
##  $ Species_full                : chr  "Trinycteris nicefori" "Lonchorina aurita" "Micronycteris hirsuta" "Phyllostomus hastatus" ...
##  $ ID                          : chr  "Naranjo" "Lonc_auri_01.08.18_1" "Ghandi" "King Lear" ...
##  $ Sex                         : chr  "M" "M" "M" "M" ...
##  $ Date                        : chr  "4/8/19" "1/8/19" "11/8/18" "2/1/19" ...
##  $ Set                         : chr  "2" "3" "2" "?" ...
##  $ Attempt                     : int  1 1 1 1 1 2 1 1 1 1 ...
##  $ Playback                    : chr  "H. frenatus" "Q. gigas" "Beetle flight" "P. pustulosus" ...
##  $ Number                      : chr  "20" "16" "17r" "2" ...
##  $ Time                        : chr  "0:00" "22:48" "23:09" "19:11" ...
##  $ Response                    : int  0 2 1 0 0 0 1 1 0 0 ...
##  $ Notes                       : chr  "0" "1" " " " " ...
##  $ Reviewed                    : chr  "" "" "" "" ...
##  $ Exclude                     : chr  "" "" "" "" ...
##  $ call_duration               : chr  "2.13" "15.735" "0.490975" "0.3785" ...
##  $ peak_freq                   : num  3788 1900 22525 858 858 ...
##  $ min_freq                    : num  1890 700 525 812 812 ...
##  $ max_freq                    : num  11148 3900 58775 2935 2935 ...
##  $ bandw                       : num  9255 3100 58200 2120 2120 ...
##  $ entropy                     : num  0.392 0.179 0.491 0.148 0.148 ...
##  $ peaks                       : num  7.25 1 7.75 1 1 1 1 1 1 1 ...
##  $ pulses_call                 : num  10.8 66.5 NA 2 2 ...
##  $ calls_callgroup             : num  1 1 NA 1 1 1.5 1 1.5 1 2 ...
##  $ sum_pulse_duration_callgroup: num  0.256 15.735 NA 0.384 0.384 ...
##  $ duration_callgroup          : num  2.13 15.735 NA 0.379 0.379 ...
##  $ call_rate_used              : num  54 15.97 NA 1.95 1.95 ...
##  $ duty_cycle                  : num  0.0394 0.9854 NA 0.1944 0.1944 ...
##  $ file_time_with_sound        : num  0.425 58.15 43.65 20.475 20.475 ...

Check and clean data

Code “call types” like “katydid” and “cicada”

#Recode a variable with the categories

g_u<- g_u %>% mutate(Call_type = if_else(Playback %in% c("A. spatulata","C. wheeleri","D. gigliotosi" ,"T. subfalcata","V. zederstedi"), "Katydid", 
      if_else(Playback %in% c( "Acla sp.","Aclodes sp.",  "Anaxipha sp."), "Cricket", 
      if_else(Playback %in% c("D. diastema", "D. ebraccatus", "H. fleishmanni", "P. pustulosus" ,"S. sila"), "Frog", 
      if_else(Playback %in% c("Q. gigas", "Z. smaragdina"), "Cicada", 
      if_else(Playback %in% c( "H. frenatus","T. rapicauda" ), "Gecko", 
      if_else(Playback %in% c("Beetle flight", "Frog hopping" , "Mouse rustles", "Moth sounds" ),"Movement cue", 
      if_else(Playback %in% c("White noise"), "White noise","Other"
      
      ))))))))


#order factor
g_u$Call_type<- factor(g_u$Call_type, levels = c( "Cricket" , "Katydid", "Cicada" , "Frog"  ,"Gecko", "Movement cue",    "White noise" ), ordered = TRUE )

unique(g_u$Call_type)
## [1] Gecko        Cicada       Movement cue Frog         Katydid     
## [6] Cricket      White noise 
## 7 Levels: Cricket < Katydid < Cicada < Frog < Gecko < ... < White noise

Note Moth sounds were only played to a few bats, should be removed for some analyses

Original response scale:

0 - no response
1 - ears twitch in time with call
2 - orients toward sound when playing
3- flies towards speaker
4 - hovers within .5 m
5- lands on speaker

Recode response with 0, 1, 4, 5, -> Make a column that codes whether bats: 0 - (did nothing (0),
1- ear twitched or came closer(1-3),
4- approached within a half meter (4),
5- or landed (5).
(collapses 1-3, these distinctions may not be important/ espcalatory across spp.)

Also make a binary column that just codes whether or not a bat landed. ->
or 0, 1, 2-3,4,5. With the idea that orient and comes closer can be conflated.

#Recode a variable with 0, 1, 4, 5 

#g_u<- g_u %>% mutate(Response_ApproachLand = if_else(Response %in% c("0","1","2" ,"3"), "0", if_else(Playback %in% c("4"), "4", if_else(Response %in%c("5"),"5", "other"))))

g_u <- g_u %>% mutate(Response_ApproachLand = case_when(
      Response == 0 ~ 0, 
      Response == 1 ~ 1,
      Response == 2 ~ 1,
      Response == 3 ~ 1,
      Response == 4 ~ 4, 
      Response == 5 ~ 5))

g_u <- g_u %>% mutate(Response_Land = case_when(
      Response == 0 ~ 0, 
      Response == 1 ~ 0,
      Response == 2 ~ 0,
      Response == 3 ~ 0,
      Response == 4 ~ 0, 
      Response == 5 ~ 1))

Palette for coloring by Call_type

# Create a custom color scale (using RColorBrewer)
library(RColorBrewer)
CallCatColors <- brewer.pal(7,"Set1")

# Palette rearranged
#CallCatColors <- c("#984EA3", "#4DAF4A", "#377EB8", "#FF7F00", "#E41A1C", "#FFFF33", "#A6A6A6") 

# Alternative palette (OG powerpoint scheme)
CallCatColors <- c("#2A769E",  "#4D6123" , "#1E846A" , "#D0892C" , "#973721", "#573AA2" , "#A3A3A3"  )

#Alternative palette 2 (modified powerpoint scheme)
#CallCatColors <- c("#FDE84C", "#4D6123", "#19657E", "#EE964B", "#8C2635", "#6D3A47", "#A6A6A6")

# View palette
pie(rep(1,7), col=CallCatColors )

names(CallCatColors) <- levels(g_u$Call_type)
colScale_CallCat <- scale_fill_manual(name = "Call type", values = CallCatColors)
colScale_color_CallCat <- scale_color_manual(name = "Call type", values = CallCatColors)
#and then add the color scale onto the plot as needed
#summarize # of individuals per species
n_species <-g_u %>% group_by(Species_full) %>%
      summarise(N_total = n_distinct(ID, na.rm=TRUE), N_Female = length(unique(ID[Sex == 'F'])), N_Male = length(unique(ID[Sex == 'M']))) %>% arrange(desc(N_total))
## `summarise()` ungrouping output (override with `.groups` argument)
n_species
## # A tibble: 15 x 4
##    Species_full              N_total N_Female N_Male
##    <chr>                       <int>    <int>  <int>
##  1 Trachops cirrhosus             13        5      8
##  2 Lophostoma silvicolum          12        2     10
##  3 Artibeus jamaicensis           10        1      9
##  4 Micronycteris microtis         10        4      5
##  5 Phyllostomus hastatus          10        0     10
##  6 Micronycteris hirsuta           7        2      5
##  7 Phyllostomus discolor           7        1      6
##  8 Gardnerycteris crenulatum       6        3      3
##  9 Lonchorina aurita               5        1      4
## 10 Lampronycteris brachiotis       3        1      1
## 11 Trinycteris nicefori            3        0      3
## 12 Lophostoma brasiliense          2        1      1
## 13 Tonatia saurophila              2        0      2
## 14 Glyphonycteris daviesi          1        0      1
## 15 Micronycteris schmidtorum       1        0      1
#sum total n bats
total_n <- sum(n_species$N_total)
total_n
## [1] 92
#write.csv(n_species, "n_species.csv")

bootstrapped means and CIs for playback for each bat spp (for faceted plots)

#  # end goal: bootsums for all bats in all playbacks
# #but blanks for n<7 in high and low columns
# #put blanks in n<2 for mean column
# 
# # n per species
# n <- g_u %>% group_by(Species_full) %>% summarize(n = length(unique(ID)))
# 
# #bootstrap playback values within each bat spp
# Boot_sums <-      g_u %>%
#       filter(Call_type !="NA") %>%
#       mutate(effect= paste(Species_full, Playback, sep="_")) %>%
# 
#  boot_ci2(d = ., y = .$Response, x = .$effect, perms = 5000, bca = F) %>%     
#       mutate(Species_Playback = effect) %>% #keep combined column in new name
#             separate(effect, into=c('Species_full', 'Playback'), sep="_", remove=F) #re separate combined column
# 
# # add n spp
# Boot_sums <- left_join(Boot_sums, n) # n per spp. 
# 
# 
# # delete bootstrapped values below minimum sample size for plotting (could also do by total bat n, not observations per each playback)
# N_min_mean <- 10 #min value for plotting 95% bats
# N_min_int <- 10 #min value for plotting boot means
# 
# Boot_sums <- Boot_sums %>%
#       mutate(low = ifelse(n.obs>=N_min_int, low,NA )) %>% 
#       mutate(high = ifelse(n.obs>=N_min_int, high, NA )) %>%
#       mutate(mean = ifelse(n.obs>=N_min_mean, mean, NA )) 
#       
#                        

Faceted bat plots

 #code for faceted plots 


facet_plot_bats <- function(d=d, ranked= FALSE, Boots = FALSE){
 
      
# sample size per species
n <- d %>% group_by(Species_full) %>% summarize(n = length(unique(ID)))

###If Boots= TRUE, calculate bootstrapped values of group
if(Boots){
      

#bootstrap playback values within each bat spp
Boot_sums <-      d %>%
      filter(Call_type !="NA") %>%
      mutate(effect= paste(Species_full, Playback, sep="_")) %>%

 boot_ci2(d = ., y = .$Response, x = .$effect, perms = 5000, bca = F) %>%     
      mutate(Species_Playback = effect) %>% #keep combined column in new name
            separate(effect, into=c('Species_full', 'Playback'), sep="_", remove=F) #re separate combined column

# add n spp
Boot_sums <- left_join(Boot_sums, n) # n per spp. 


# delete bootstrapped values below minimum sample size for plotting (could also do by total bat n, not observations per each playback)
N_min_mean <- 10 #min value for plotting boot means
N_min_int <- 10 #min value for plotting bootstrapped 95% intervals

Boot_sums <- Boot_sums %>%
      mutate(low = ifelse(n.obs>=N_min_int, low,NA )) %>% 
      mutate(high = ifelse(n.obs>=N_min_int, high, NA )) %>%
      mutate(mean = ifelse(n.obs>=N_min_mean, mean, NA )) 
}            
      
      
           
#Set factor levels for species from high to low mean response for plotting (doesn't apparently work)   
Species_levels <- levels(fct_reorder( d$Species_full, -d$Response, .fun='mean'))
      
# various failed attempts to make the facet_labels read: "italic(species name ) n= Sample size"
      #part of the problem is that I haven't been able to reorder the facets (by changing the factor levels) outside the facet_wrap() command
            #appender <- function(string, suffix) paste0(string, suffix)
            #appender <- function(string) paste0(string, distinct(subset(Boot_sums, Boot_sums$Species_full == string), n))
            # returns all n, not the n of facet
            #p + facet_wrap(~am, labeller = as_labeller(appender)) 
            #label facets
             lev <- n$Species_full
            lab <-      paste(n$Species_full, " (n = ", n$n , ")", sep = ""   )
            names(lab) <- lev

            
#make dotplot for bats       
         d %>%
           filter(Call_type !="NA") %>%
           mutate(Species_full = factor( Species_full, levels = Species_levels)) -> d
    
#Why don't these reorder facets? infuriating
     d$Species_full <- fct_reorder( d$Species_full, -d$Response, .fun='mean')

     
                  
#reorders the x-axis from high to low mean response if ranked = TRUE
 if(ranked){
    plot <-   
          d %>%
           mutate(Playback = fct_reorder(Playback, Response, .fun='mean')) %>%
          ggplot( aes(x = reorder(Playback, -Response), y = Response))  
    }else{
                plot <-   d %>%
                      ggplot(aes(y=Response, x = Playback))
    }
     
     plot <- plot +
           facet_wrap(~factor(Species_full,
                              levels = Species_levels), #order facets by mean response
                      #facets = ~Species_full,
                      ncol = 1,
                      #labeller = as_labeller(appender),
                      #labeller = as_labeller(lev),
                      scales = ifelse(ranked, "free_x", "fixed")) +
           geom_dotplot( aes(color = Call_type,
                              fill = Call_type ) ,
                         binaxis = "y",
                         stackdir = "center",
                         method = "histodot" , # fixed bin width dotdensity is default
                         binwidth =1/40, #max bin width - proportion of range of data
                         alpha = .7,
                         stackratio = 0.5, # doesn't  change dot size #.5 .9 for seperate good for together plot
                          dotsize = 6 ,  #(previously not specified) #2 good for collective plot 8 or 12 good for single spp. plots
                         # width = .2
                         # right = TRUE
                        ) +
                        scale_fill_manual( 'Call type', values = CallCatColors) +
                        scale_color_manual('Call type', values = CallCatColors) +
                        # colScale_CallCat + # fill palette
                        # colScale_color_CallCat + #color palette
                       
            coord_cartesian(ylim = c(0, 5)) + #standardize y axis
          
          theme(panel.grid.major = element_blank(),
                panel.grid.minor =  element_blank(),
                  panel.background = element_rect(fill = 'white'),
                  axis.line = element_line(colour = "black"),
                  text = element_text(color = "black", size = 10),
                  plot.title = element_text(face = "italic", size = 11),
                  legend.position = "top",
                  axis.title.x = element_blank(),
                  axis.text.x = element_text(angle = 63, hjust=1, face = "italic"),
                  axis.text = element_text( color = "black", size = 10), 
                  #axis.title = element_text(size = 16),  
                  )  + 
          #dark_theme_noline +
          xlab("Playback") +
          ylab("Response scores") 
           
     #add bootstrapped means if Boots = TRUE
     if(Boots){
     plot = plot +    geom_point(data = Boot_sums,
                                          aes( x= Playback,
                                               y = mean),
                                          color = "black",
                                          shape = 4,
                                          size = 1 ) +
           
           
      #if more than n individuals, add bootstrapped 95% CI's to plot
         geom_errorbar( data = Boot_sums,
                                aes( x=Playback,
                                     ymin = low,
                                     ymax =high,
                                     y = NULL ),
                                color = "black",
                                width=.2,
                                size=.25)
     }
     
      # #make pdfs of plots (85, 129, or 177 mm width to fit in 1-3 columns dimensions)  
       # ggsave("testplot.jpg",  width = 177, height =800, units = "mm")
                 plot = plot +   ggsave(filename = paste("faceted_", ifelse(ranked,"ranked_", "unranked_"), ifelse(Boots,"CIs_", "no_CIs_"), "all_bats_", "responses_to_playbacks", ".jpg", sep = ""), width = 177, height =800, units = "mm" )
                 
      #add plot to list           
     list(plot)
      
}

unranked_plot <- facet_plot_bats(d = g_u, ranked = FALSE, Boots = FALSE)
## `summarise()` ungrouping output (override with `.groups` argument)
unranked_plot
## [[1]]

# 
# ranked_plot <- plot_bats(d= g_u, ranked = TRUE)
# ranked_plot


#TO DO: 
#make facet titles italic and add N's (use hum_names strategy below?)
#collapse response scale 3 into 2 (not apparently informative)
#Low priority: ranked isn't working right (uses same order for all facets). Fix? 
#pull out dotplot variables into function so can change for each plot as makes sense
#low priority: put Bootsums into function so not hanging out
#put the other plotting variables into the function


#respace dots

####### IGNORE ME. more workspace for figuring out how to add sample sizes to label facets

# #Create an object using as_labeller(). If the column names begin with a number or contain spaces or special characters, don't forget to use back tick marks:
# # Necessary to put RH% into the facet labels
# hum_names <- as_labeller(
#      c(`50` = "RH% 50", `60` = "RH% 60",`70` = "RH% 70", 
#        `80` = "RH% 80",`90` = "RH% 90", `100` = "RH% 100"))
# #Add to the ggplot:
#     ggplot(dataframe, aes(x = Temperature.C, y = fit)) + 
#         geom_line() + 
#         facet_wrap(~Humidity.RH., nrow = 2, labeller = facet_names)

#make list that reads "Species full = "Species full (n= N)

# facet_names = n %>%
#       mutate(facet_names = paste(Species_full, " = " , Species_full, " (n = ", n, ")" ,  sep = ""  )) %>% select(facet_names)
# 
# 
# # figure out how to add quotes around Species full (n= n_ )
#   n<-      n %>%
#        mutate(facet_names = paste(Species_full, " (n = ", n , ")", sep = ""   )) =
# 
# # labeller =  as_labeller(facet_names)
# labeller = label_bquote(italic(.(sex))~subjects))
# 
# labeller = label_bquote(italic(.(Species_full))~subjects))
# # .
# length(unique(.$ID)
# abeller = labeller(cyl = 
#     c("4" = "condition: 4",
#       "6" = "condition: 6",
#       "8" = "condition: 8"))       

  #labeller = label_bquote(italic(.(factor(Species_full, levels = Species_levels)))),
                      #labeller =  as_labeller(facet_names),

# lev <- n$Species_full
# lab <-      paste(n$Species_full, " (n = ", n$n , ")", sep = ""   )
# 
# names(lab) <- lev
#      
# 
# label_facet <- function(original_var, custom_name){
#   lev <- levels(as.factor(original_var))
#   lab <- paste0(lev, ": ", custom_name)
#   names(lab) <- lev
#   return(lab)  
# }
# 
# facet_wrap(~cyl, labeller = labeller(Species_full = label_facet(original_var = n$Species_full, custom_name = "grouping")))
# 
# facet_wrap(~cyl, labeller = labeller(cyl = label_facet(df$cyl, "grouping")))
###### IGNORE ABOVE ################

plot with only top

Individual plots, one bat per plot, ranked and unranked x axis

#Y= response, X= playback, within bat spp. boot_ci()

# store dataframes and plots in lists to be able to make a multiplot or grid object- see how plots look together 
# remove all the legends
#try making a multiplot situation (ggarrange from egg package, ?align_plots from the cowplot package?)
#make a version where X isn't reordered every time, so just one X on bottom,
#

#response score- 3 has so few data points- doesn't taper with more- combine 3 and 2?

#clean this up (get rid of excess lists and whatnot)





#pull unique species names
Spp <- unique(g_u$Species)

#initialize lists for plots and boot output


####fixes#####
# pull out ranked from function, test if function works on its own x
#try adding in lists, internally x
#check what output of mutate is
# maybe pull out x and y? response and playback? to make function more universal?
##############

#make seperate plots and bootstrapped 95% cis for each bat
# ranked = X axis ranked from greatest to smallest mean response
# N_min_mean = minimum # for plotting boot means
# N_min_int = minimum # for plotting boot intervals
plot_bats_separate <- function(d=d, ranked= TRUE, N_min_mean=10, N_min_int=10) {
      
      #initialize lists for plots and boot output
      myplots <- list()
      bootsummary <- list()
      
        #loop through bat spp
        for (i in 1:length(Spp)) {
      
  #pull out one bat species at a time          
  data <-      d %>%
                filter(Call_type !="NA") %>%
                filter(Species ==Spp[i])
  
   # # of bats in group 
   n <- length(unique(data$ID))
  

    #get bootstrapped means and 95% confidence intervals for the response to each playback in n> 10

   if(n >=N_min_mean){
           Boot_sums <- boot_ci2(d = data, y = data$Response, x = data$Playback, perms = 5000, bca = F) %>%
                 rename(Playback = effect)

   }
          
           
      #store bootsums in list     
          #  bootsummary[[i]] <- Boot_sums
      
      #make dotplot for bats       
     plot <-     data %>%
                 filter(Call_type !="NA") 
  
        
#reorders the x-axis from high to low mean response if ranked = TRUE
 if(ranked){
    plot <-   plot %>% mutate(Playback = fct_reorder(Playback, Response, .fun='mean')) %>%
          ggplot( aes(x = reorder(Playback, -Response), y = Response))  
    }else{
                plot <-   plot %>%
                      ggplot(aes(y=Response, x = Playback))
    }
     
     plot <- plot +
           geom_dotplot( aes(color = Call_type,
                              fill = Call_type ) ,
                         binaxis = "y",
                         stackdir = "center",
                         method = "histodot" , # fixed bin width dotdensity is default
                         binwidth =1/40, #max bin width - proportion of range of data
                         alpha = .7,
                         stackratio = 0.9, # doesn't  change dot size #.5 good for together plot
                          dotsize = 8 ,  #(previously not specified) #2 good for collective plot 12 good for single spp. plots
                         # width = .2
                         # right = TRUE
                        ) +
                        colScale_CallCat + # fill palette
                        colScale_color_CallCat + #color palette
              
            coord_cartesian(ylim = c(0, 5)) + #standardize y axis
          
          theme(panel.grid.major = element_blank(),
                panel.grid.minor =  element_blank(),
                  panel.background = element_rect(fill = 'white'),
                  axis.line = element_line(colour = "black"),
                  text = element_text(color = "black", size = 10),
                  plot.title = element_text(face = "italic", size = 11),
                  legend.position = "right",
                  axis.title.x = element_blank(),
                  axis.text.x = element_text(angle = 63, hjust=1, face = "italic"),
                  axis.text = element_text( color = "black", size = 10), 
                  #axis.title = element_text(size = 16),  
                  )  + 
          #dark_theme_noline +
          xlab("Playback") +
          ylab("Response scores") +
          ggtitle(paste(unique(data$Species_full)," (n = ", length(unique(data$ID)), ")", sep ="") )   #makes title current group, n is number of bats      
        
           if(n >=N_min_mean){
                 plot = plot + geom_point(data = Boot_sums,
                                          aes( x= Playback,
                                               y = mean),
                                          color = "black",
                                          shape = 4,
                                          size = 1 )
           }
           
      #if more than n individuals, add bootstrapped 95% CI's to plot
          if(n >=N_min_int){
                plot = plot + geom_errorbar( data = Boot_sums,
                                aes( x=Playback,
                                     ymin = low,
                                     ymax =high,
                                     y = NULL ),
                                color = "black",
                                width=.2,
                                size=.25)
          }
          
      # #make pdfs of plots (85, 129, or 177 mm width to fit in 1-3 columns dimensions)            
                 plot = plot +   ggsave( filename = paste("individual_", ifelse(ranked,"ranked_", "unranked_"), unique(data$Species_full), "responses_to_playbacks", ".jpg", sep = ""), width = 177, height = 55, units = "mm" )
                 
      #add plot to list           
      myplots[[i]] <- plot
      
    
      #print(plot)
     # print(Boot_sums)
        }
names(myplots) <- Spp
myplots   
}


#unranked plots
unranked_individ_batplots <- plot_bats_separate(d = g_u, N_min_mean = 10, ranked = FALSE)
#unranked_individ_batplots

#ranked plots
ranked_individ_batplots <- plot_bats_separate(d = g_u, N_min_mean = 10, ranked = TRUE)
ranked_individ_batplots
## $Trin_nice

## 
## $Lonc_auri

## 
## $Micr_hirs

## 
## $Phyl_hast

## 
## $Loph_silv

## 
## $Glyp_davi

## 
## $Gard_cren

## 
## $Trac_cirr

## 
## $Micr_micr

## 
## $Arti_jama

## 
## $Tona_saur

## 
## $Phyl_disc

## 
## $Loph_bras

## 
## $Lamp_brach

## 
## $Micr_schm

combined plot with ranked x axes

#make plot with aa bats and ranked X axis 
ggpubr::ggarrange(plotlist = ranked_individ_batplots[1:15], ncol = 1, common.legend = TRUE)

ggsave("AA_test_plot_allbats_long.jpg", width = 177, units = "mm", height = 800)

Add plots with 1 plot per playback (~line 950)